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We improve the calculation of the equation of state for two flavor QCD by simulating on A^t = 6 lattices at 
appropriate values of the couplings for the deconfinement/chiral symmetry restoration crossover. For am, = 
0.0125 the energy density rises rapidly to approximately 1 GeV/fm'' just after the crossover(m^/mp ~ 0.4 at this 
point). Comparing with our previous result for A^'t = 4 |lj, we find large finite Nt corrections as expected from 
free field theory on finite lattices. We also provide formulae for extracting the speed of sound from the measured 
quantities. 



1. INTRODUCTION 

In order to show the existence of the Quark- 
Gluon Plasma (QGP) in the aftermath of upcom- 
ing heavy-ion collision experiments at RHIC and 
CERN or to understand the dynamics of the QGP 
in the early universe, one needs as input, among 
other things, the equation of state for QCD, i.e. 
the energy density and the pressure as a function 
of temperature and quark mass. 

We are continuing our program of computing 
the equation of state for two flavor QCD. Last 
year we reported results for Nt = 4 [0 , and we are 
now working at A^t = 6. A similar program for the 
pure gauge theory is being pursued by the Biele- 
feld group 1^. The Nt — 6 simulations represent a 
significant increase in computational cost due to 
the increased lattice size, smaller quark masses, 
and corresponding smaller simulation step sizes. 
Step size (At) errors induced by the approximate 
integration of the gauge field equations of mo- 
tion are much more significant at A'^t = 6 than at 



Nt — 4. They are handled by extrapolation of 
observables to At = 0, and by running at small 
step sizes, with a corresponding large increase in 
the cost of simulation. 

We have surveyed the gauge coupling and 
quark mass plane for two flavor QCD in a region 
relevant to the nonzero temperature crossover in 
order to measure the nonperturbative pressure by 
integration. The interaction measure is also cal- 
culated and together with the pressure it yields 
the energy density. 

2. THEORY 

A Euclidean N^ x Nt lattice with periodic 
boundary conditions has a temperature T and 
volume V 

V = N!a^ 

l/T = Nta, (1) 

where a is the lattice spacing. Thermodynamic 
variables are derivatives of the partition function 
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Z . In particular, the pressure p and energy den- 
sity £ are given by 



p 91og Z 
T " dV 



and 



eV = - 



dlogZ 

Wm 



(2) 



(3) 



The methods for computing the pressure and 
integration measure are discussed in Ref . jlj . The 
pressure is found by integrating either the plaque- 
tte or 
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The interaction measure is given by 

rV_ _ 1 d\ogZ _ d\ogZ 
'Y ~ ~Td{l/T) ~ dV 

oat Otts mog a 

where as and at are the spatial and temporal lat- 
tice spacings. The scale dependence in the case 
of QCD with dynamical quarks leads to 
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The derivatives are the usual (3 function and the 
anomalous dimension of the quark mass. In the 
above the subscript "sym" refers to symmetric 
lattices with Nt — which are used for vacuum 
subtraction. 

The knowledge of the non-perturbative pres- 
sure and interaction measure allows us to com- 
pute other bulk quantities. The energy and en- 
tropy s become 



e ^ I + 3p 
,sT = I + 4p . 



(8) 



3. SOUND SPEED 

Acoustic perturbations travel in the system 
with a speed Cs'. 

- = - (9) 

C2 dp 

One has to take the derivative keeping the phys- 
ical quark mass fixed, i.e. on the line of con- 
stant physics. Unfortunately, we know best only 
the variations of the energy density and pressure 
along lines of constant bare parameters. In or- 
der to measure the correct sound speed one has 
to use the P function to map the changes in bare 
parameters to physical changes of temperature. 

For QCD with zero chemical potential the only 
varying quantity is the temperature T. In the 
real world the masses of the quarks do not change 
and derivatives are taken at fixed quark mass. 
The volume is infinite and it is divided out of the 
equations that consider densities. Hence 

1 drl^.^Tr/TOp 



dp I 

dT \V,ra^jrap 



(10) 



Henceforth we drop the |m„/mp reminder. With 
only T varying, the fundamental relation of ther- 
modynamics becomes 



/(T) = e{t) - r,s(r) = -p 



(11) 



or 



(12) 



e[T)^Ts{T)^p{T). 

Then, 
dp 

— = Ts\T) + s{T) p'[T) = Ts'iT), (13) 

where we have utilized a Maxwell relation for the 
entropy s: 



(14) 



In the Maxwell relation we interchanged pressure 
and free energy according to 



PV , „ -fv 

— log Z = ^ . 
Hence, 

1 _ Ts'{T) _ rflog(s) 
c| ~ p'{T) ^ dlogT 



(15) 



(16) 
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where we have used Eq. (|lj) again. 

However, lattice simulations are always done at 
finite volume, and the volume varies with the lat- 
tice spacing a. To get the correct value one must 
take the derivative with respect to temperature 
with the volume constant. 

The derivative with respect to temperature is 
problematic. It requires asymmetric lattice spac- 
ing and leads to expressions with asymmetry co- 
efficients. These in turn, are poorly known in 
the regime of bare couplings, where simulations 
are currently feasible. It would be advantageous 
to find an expression that does not involve the 
asymmetry coefficients. 

Using Eq. (||), one can give the sound speed 
161) with the interaction measure: 



1 



1 dl 



On the other hand, Eq. 

r T aiog z 1 



V 91og a J 



dT 



(17) 

allows us to write 

(18) 



Now, the T derivative is taken at constant V . 
Therefore, 



ci s 



laiogz T d[^] -' 

V dloga V dT 
The last derivative can be given as 

dT T dloga 



(19) 



(20) 



which can be seen by expanding the derivatives 
with respect to a and T as derivatives with re- 
spect to at and as: 
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(21) 



We also used the definition of energy density, Eq. 
(||), in the form 



eV _ dlogZ 



(22) 



Then Eq. ( [19| ) becomes 
T d[ 



1 



1 



T i 



V dlog a 



(23) 



Using sT = e + p and / = £ — 3p this can be 
expressed in a rather compact form 
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4e- 



T d[^] 
V dlog a 



(24) 



If the system is conformally invariant the energy 
density does not depend on the scale and the 
derivative term in (p4r) disappears: 
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= 3, 



P 



(25) 



the relativistic result for a free gas. 

For QCD with dynamical quarks an explicit 
form is 



1 



1 



ea^ 



■ pa^ 



dea'^ 



dloga 'd(6/g^) 
^dm„a , dea^ 



(26) 
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The derivatives have to be taken on a line of con- 
stant physics. 

4. SIMULATIONS 

We have measured (□) and {ipip') on asymmet- 
ric lattices with Nt = 6 and Ng = 12, and on 
symmetric lattices with Nt — Ng — 12. 

The couplings in the simulations are appropri- 
ate for the Nt = 6 crossover |^. At aniq = 0.0125 
and 0.025 we varied between 5.37 and 5.53 
and 5.39 and 5.53, respectively. These runs will 
allow us to extrapolate results to zero quark mass. 
At = 5.45 and 5.53, we have varied anig be- 
tween 0.01 and 0.1 and 0.0125 and 0.2, respec- 
tively. These runs also allow us to extrapolate to 
zero quark mass, but in addition they provide a 
cross check on the integrations over and give 
information on the equation of state over a wide 
range of quark masses. 

Since we are simulating with two flavors of 
Kogut-Susskind fermions, the simulations are 
performed with the refreshed molecular dynam- 
ics R algorithm H. For the Nt = 6 lattices we 
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ran at least 1800 trajectories after 200 trajecto- 
ries for thermalization. On the symmetric lattices 
we performed at least 800 trajectories after 200 
trajectories for thermalization. Each trajectory 
had unit length in simulation time. 

The R algorithm induces an error in the ob- 
servables with a leading term proportional to the 
square of the step size, in simulation time, used 
to integrate the equations of motion of the gauge 
fields. In practice, the error is small for each ob- 
servable. However, the errors are different on the 
symmetric and asymmetric lattices, so they do 
not cancel from the vacuum subtractions. More- 
over, in many instances the errors are the same or- 
der of magnitude as the vacuum subtracted quan- 
tities. Therefore, they must be eliminated. 

In most cases, the step size errors are elimi- 
nated by extrapolation. For each gauge coupling 
and quark mass several simulations are run with 
two or more step sizes. If the step sizes are small 
enough, observables depend linearly on the step 
size squared. Once this linear dependence is ob- 
served, the quantities are extrapolated to zero 
step size. As a rule of thumb, the step size in 
the R algorithm should be less than or approx- 
imately equal to aniq. Roughly speaking, this 
is because in the integration of the gauge mo- 
mentum, the fermion "force" to lowest order is 
proportional to l/amq. Thus the step taken in 
simulation time to update the momentum should 
be ~ arriq (or smaller) to keep the change in the 
momentum less than 0(1). Thus we have used 
0.007 < At < 0.015 and 0.015 < At < 0.03 
for the runs with amq — 0.0125 and 0.025, re- 
spectively. For the larger quark mass runs at 
6/52 = 5.45 and 5.53, we used 0.02 < At < 0.03. 
Finally, for the smallest quark mass simulation, 
aruq = 0.01 (6/.g2 = 5.45), we use At = 0.005. 

5. EXTRAPOLATIONS 

In Fig. 1^ we show the plaquette dependence 
on step size squared for aniq = 0.0125. Similar 
results hold for arUq = 0.025. Generally, the ef- 
fects are larger on the symmetric lattices and at 
smaller For example, at = 5.39, lin- 

ear behavior sets in for much smaller step size on 
the cold lattice than on the hot lattice. In fact. 



1.70 



1.65 



A 
□ 
V 



1.60 



1.55 

0.0000 0.0001 0.0002 

Figure 1. The plaquette as a function of At^ for 
aniq = 0.0125. Lines are to guide the eye and are 
not fits. Labels for the couplings refer to the hot 
lattices (octagons). 



for At"^ > 0.0001, one would conclude that the 
system is in the confined phase; only at smaller 
step size is a clear separation visible. From Fig. |l|, 
the following is a reasonable extrapolation proce- 
dure. First, on the symmetric lattices for 5.39 < 
6/5^ < 5.43, use only the smallest two step sizes 
to extrapolate to At'^ = 0. At 6/5^ = 5.37 we 
effectively assume the system is in the cold phase 
and take the values at At = 0.007 as the zero 
step size extrapolations. For > 5.43 we use 
all At values to do the extrapolations. On the hot 
lattices we do the following. For each where 
there are multiple step sizes, we use all of them 
to extrapolate to At'^ = 0. For 6/3^ = 5.40 and 
5.42 we interpolate the slope from the neighbor- 
ing points and use that along with the point at 
At — 0.01 to extrapolate to zero step size. Note, 
for 6/g^ > 5.43, the slopes are essentially zero. 
Therefore, for 6/g^ > 5.43 where only one mea- 
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surement is available, we take that value as the 
zero step size value. 

For (^jjil>) the situation is similar to the pla- 
quette with the following exceptions. First, after 
vacuum subtraction, the relative step size errors 
are not as significant. Second, on the hot lattices 
for 6/(7^ > 5.39, we find the slopes with respect 
to Ai^ are zero within errors, so we take the zero 
step size {i^i^) to be the value measured at the 
smallest At for each of these couplings. 

For aruq = 0.025, the situation is similar, but 
the smallest step size runs are still incomplete as 
of this writing. 

Because the cold lattices vary smoothly with 
6/5^, we made cold runs at values of 6/(7^ sepa- 
rated by 0.02, while the hot runs were separated 
by AG/g'^ — 0.01 near the transition. Cold ob- 
servables at the other couplings were obtained 
by interpolation. For 6/(7^ < 5.47, quadratic 
fits to the zero step size plaquette and had 
= 2.69 and 3.18 with three degrees of freedom 
respectively. These fits were used for the interpo- 
lated values. To get the symmetric observables at 
6/(7^ = 5.53, one can either extrapolate in 6/(7^ 
for fixed am,q, or extrapolate in aniq for fixed 
6/g2. Both give the same result within errors, 
and we simply take the value from extrapolation 
in arUq as it has a much smaller error. 

6. PRESSURE 

The results for (t/'V') ^^J^d (□) from the previous 
section can now be integrated to yield the pres- 
sure. 

To begin consider (V"/') &s a function of aniq. 
Using Eq. (||), we find the pressure as a func- 
tion of aniq (at 6/5^ = 5.45 and 5.53). The 
result is shown in Fig. ^. We also want the 
pressure at arUq = which is found by setting 
ipipiO) = and continuing the integration to 
aniq — 0. At — 5.45, a linear fit to the 

data for aniq < 0.025 that is constrained to go 
through the origin has = 4.9 for three degrees 
of freedom. For = 5.53, we only have mea- 
surements at aniq = 0.0125, 0.025, and 0.05 for 
which the linear fit does not work well. However, 
a quadratic fit constrained to go through the ori- 
gin has — 0.49 for one degree of freedom. The 



pressure at aniq — calculated in this manner is 
also shown in Fig. ||. 
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Figure 2. The pressure from integration of 
< ipip > with respect to aniq. The bursts are ex- 
trapolations to zero quark mass. 



Next we integrate the plaquette with respect 
to 6/(7^ to obtain the pressure as a function of 
6/(7^ at fixed aniq. The result for aniq — 0.0125 
is shown in Fig. The pressure rises smoothly 
through the crossover region as it must. The re- 
sults from the quark mass integrations are also 
shown in Fig. ^, and are in agreement with the 
integration. This is a good check on our 
analysis, since for the most part the two integra- 
tions are independent. In particular, the integra- 
tion with respect to is sensitive to the step 
size extrapolations while the amq integration is 
not. 



7. INTERACTION MEASURE 

The interaction measure is given by Eq. (^. 
For the f3 function we use our previous result cal- 
culated from the tt and p masses at various values 
of 6/g^ and aniq The result is shown in Fig. IJ. 
It rises sharply from zero in the cold phase to 
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Figure 3. The pressure from integration of (□) 
with respect to {anig — 0.0125). The crosses 
are from the {ip'4') integration. 
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Figure 4. The interaction measure for arriq 
0.0125. 



some maximum and then begins to drop off. For 
a noninteracting plasma, / is zero since s — 3p 
which is certainly not the case at the largest value 
of 6/(7^ in our simulations. From asymptotic free- 
dom we expect the system to approach a nonin- 
teracting plasma at very high temperature. 

The interaction measure at aruq — depends 
only on the plaquette since the anomalous dimen- 
sion of the quark mass is zero at arUq = 0. To 
extrapolate the plaquette to zero quark mass we 
use the fact that its slope is just its correlation 
with ipip, 



d{0) 
d{amq 



(27) 



Since rpip is discontinuous at the origin in the bro- 
ken phase and continuous in the chirally symmet- 
ric phase, we expect a cusp at the origin for the 
cold lattices and zero slope for the hot lattices (l). 
At = 5.45 a linear fit on the cold lattices for 
aniq < 0.05 gives = 1-6 with two degrees of 
freedom. On the hot lattice a quadratic fit con- 
strained to zero slope at the origin gives = 0.74 
with two degrees of freedom. At = 5.53 the 
situation is less satisfactory. For the cold lattices. 



a quadratic fit to the data with aruq < 0.1 has 
— 0.09. Note the smallest quark mass here 
is 0.025, and we also use this fit for the cold ob- 
servables at aniq = 0.0125. On the hot lattices 
the data seem to reach a maximum at nonzero 
quark mass, so we take the measurement at the 
smallest quark mass to be the extrapolated value. 
Similar behavior at G/g^ = 5.53 was observed in 
our earlier Nt — 4: study, which may indicate a 
systematic error. Naively one expects finite vol- 
ume effects to order the lattice which is opposite 
to the observed behavior. The data and the fit 
results are shown in Fig. ||. 

8. EQUATION OF STATE 

In Fig. ^ we show the Nf = 6 equation of state 
as a function of 6/17^ for aniq = 0.0125. The zero 
quark mass extrapolations are also shown. At 
the largest value of 6/17^, e is still much larger 
than 3p; for aniq = 0, 3p is approximately 75% 
of e. On the other hand, after a rapid rise, e/T'^ 
more or less levels off for 6/(7^ > 5.43. This corre- 
sponds to an energy density of about 1 GeV/fm'^ 
at 6/g^ — 5.43 (we use the p mass to convert 6/g^ 
to temperature). mT^/mp at this point is slightly 
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Figure 5. The plaquette as a function of the 
quark mass. The lower (upper) curves are for 
— 5.45(5.53). The bursts are results from 
fits to the data except for the hot curve (oc- 
tagons) at = 5.53 where the measured value 
at aruq = 0.0125 is taken as the zero quark mass 
result. 



less than 0.4. 

In Fig. 1^ we compare the Nt = 4 and 6 equa- 
tions of state. There is a large finite size ef- 
fect as expected from the free field results on 
finite lattices (also shown). For example, p/T'^ 
for aniq — differs by about 15%. There ap- 
pears to be a sizable quark mass effect as well 
(the Nt — 4: results are for amq — 0.1 and 0.025). 
For Nt — 6 the approach to the Stefan-Boltzmann 
law is unclear since we do not have data at very 
high temperatures. The prominent peak in e/T'^ 
for iVt = 4 just after the crossover is much smaller 
at Nt — 6. We have plotted the results versus 
temperature by using the p mass to set the scale. 
The crossover temperature is around 150 MeV 
and is rather insensitive to Nt and aruq. In a re- 
cent preprint, Asakawa and Hatsuda have pointed 
out that many features of this equation of state 
are constrained by fundamental thermodynamic 
relations [pi. 



Figure 6. The Nt — 6 equation of state for 
amq = 0.0125 and extrapolations to arriq ~ 
O(bursts). The upper curve is e/T'^. 



9. SOUND SPEED 

The sound speed is a quantity that depends on 
the second derivative of the partition function. 
Therefore it is more difficult to get its value than 
the values for the thermodynamic variables dis- 
cussed so far. Even if our formulation can avoid 
the asymmetry coefficients there is an awkward 
derivative of the energy density in the formula. 

To measure the change in the mass as accu- 
rately as possible we performed a set of simula- 
tions at amq — 0.09 with Nt = 4 lattices in addi- 
tion to our old data at amq = 0.1. The result is 
shown in Fig. ^. Close to the transition the error 
in the derivative of the energy density overwhelms 
our data and we have large error bars. For large 
temperatures, we see that the speed approaches 
the ideal gas value. 

10. CONCLUSION 

We have calculated the equation of state for 
QCD with two flavors of quarks on iVt = 6 lat- 
tices. The algorithm used to generate gauge con- 
figurations introduces a step size error in observ- 
ables that must be removed by extrapolation. 
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Figure 7. Comparison of the equation of state 
for TVj = 4 (solid lines) and 6 (dashed lines). The 
results shown are for aniq = 0.0125 (diamonds), 
0.025 (octagons), and 0.1 (squares). Bursts are 
extrapolations to aruq = 0. The horizontal lines 
give the Stefan-Boltzmann law for Nt — 4, 6, and 
the continuum (lowest line). 



The added computational cost is significant. For 
arUq = 0.0125, we find the energy density just af- 
ter the crossover to be roughly 1 GeV/fm'^, and 
for T almost twice the critical value, the high- 
est that we simulated, three times the pressure is 
only 60% of the energy density. We find large ef- 
fects of nonzero lattice spacing from comparison 
with results at Nt = 4, as expected from free field 
theory. 

After the completion of runs at am^ = 0.025, 
we will complete our extrapolation to zero quark 
mass. It still remains to eliminate remaining lat- 
tice size effects, include the effects of the strange 
quark, eliminate any effects of finite volume, and 
include the effect of nonzero net quark density. 

This work was supported by the US DOE and 
NSF. Computations were done at the San Diego 
Supercomputer Center, the Cornell Theory Cen- 
ter, and Indiana University. 
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Figure 8. The inverse sound speed squared for 
the Nt — 4: system at niqa = 0.1. 
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